%% Compare only DMS test period data (Sensory), with Mental rotation data
%  during the response period, which should be pure mental responses.
%
% [mfr_sensory mfr_imagery]=generate_mfr_fmri_analysis( seed)

movement_trial_start_times1 = [ 490000 535000  580000 625000 670000 715000 760000 805000];
movement_trial_start_times2 = movement_trial_start_times1 + 415000;

seed = 0;
[mfr_sensory mfr_imagery]=generate_mfr_fmri_analysis2( seed, movement_trial_start_times1, movement_trial_start_times2);

seed = 100;
[mfr_sensory(:,2) mfr_imagery(:,2)]=generate_mfr_fmri_analysis2( seed, movement_trial_start_times1, movement_trial_start_times2);

seed = 200;
[mfr_sensory(:,3) mfr_imagery(:,3)]=generate_mfr_fmri_analysis2( seed, movement_trial_start_times1, movement_trial_start_times2);

seed = 300;
[mfr_sensory(:,4) mfr_imagery(:,4)]=generate_mfr_fmri_analysis2( seed, movement_trial_start_times1, movement_trial_start_times2);

seed = 400;
[mfr_sensory(:,5) mfr_imagery(:,5)]=generate_mfr_fmri_analysis2( seed, movement_trial_start_times1, movement_trial_start_times2);


%% Analyze data area by area
% s_gt_im = zeros(size(mfr_sensory,1),1);
% im_gt_s = s_gt_im;

% mfr of area different on average across subjects?
clear p h
i=1;
[first_neuron(i) last_neuron(i)]=get_group_indices('V1','TCs')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Visual Input Thal';

i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('V1','p23')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Visual Features E2/3';


i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('V1','p6(L4)')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Visual Features E5/6';


i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('P','p23')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Sequence E2/3';


i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('P','p6(L4)')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Sequence E5/6';


i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('WM','p23')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Working Memory E2/3';


i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('MTS','p23M')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Decision E2/3';


i=i+1;
[first_neuron(i) last_neuron(i)]=get_group_indices('M1','p5(L56)')
mfr_sensory_by_area(i,:)=mean(mfr_sensory(first_neuron(i):last_neuron(i),:));
mfr_imagery_by_area(i,:)=mean(mfr_imagery(first_neuron(i):last_neuron(i),:));
[p(i) h(i)]=ranksum(mfr_sensory_by_area(i,:), mfr_imagery_by_area(i,:));
name{i} = 'Motor E5/6';

mean_diff = [mean(mfr_imagery_by_area,2)- mean(mfr_sensory_by_area,2)];

sem = nansem((mfr_imagery_by_area-mfr_sensory_by_area)');

bar(mean_diff,'w','LineWidth',3);
hold on;
errorbar(mean_diff,sem,'k.','LineWidth',2)
drawnow
set(gca,'XtickLabel',name,'FontSize',14)
xlabel('Neural population','FontSize',16);
ylabel('Mean firing rate difference (Hz)','FontSize',16);
rotateXLabels( gca, 45 ); % From MATLAB central file exchange.
set(gca,'FontSize',12)

axis([.5 8.5 -1 3.0 ])

j=1;
comparisons = [];
ps = [];
for i = 1:length(name)
    if p(i) < .05;
        comparisons{j} = [i i-.15];%[i-.15 i+.15];
        ps(j)=p(i);
        j=j+1;
    end
end
%ps(ps<.01)=.011;  % Only want single asterisks; 
sigstar(comparisons,ps)


%% Analyze firing cell overlap area by area.
disp('% neurons active');
disp('% of ACTIVE active only in imagery')
disp('% of ACTIVE active only in sensory')
disp('% overlap');
disp('r of imagery and sensory population firing rates; ACTIVE neurons')

Nsubjects = 5;
for i = 1:length(name)

    disp(name{i});
    for s = 1:Nsubjects

        active_imagery = mfr_imagery(first_neuron(i):last_neuron(i),s)>1;
        active_sensory = mfr_sensory(first_neuron(i):last_neuron(i),s)>1;

        N_neurons = length(active_imagery);

        active = active_sensory|active_imagery;
        N_active = sum(active);

        pct_active(s) = N_active/N_neurons*100;

%         disp('% of ACTIVE neurons active in imagery');
%         sum(active_imagery)/N_active*100
%         disp('% of ACTIVE neurons active in sensory');
%         sum(active_sensory)/N_active*100


        pct_imagery(s) = sum(active_imagery&~active_sensory)/N_active*100;

        pct_sensory(s) = sum(active_sensory&~active_imagery)/N_active*100;

        overlap = active_sensory&active_imagery; 
        pct_overlap(s) = sum( overlap )/N_active*100;

        mfr_area_imagery = mfr_imagery(first_neuron(i):last_neuron(i),s);
        mfr_area_sensory = mfr_sensory(first_neuron(i):last_neuron(i),s);
        r(s) =corr(mfr_area_imagery(active),mfr_area_sensory(active));

        format compact 
        format short
        disp([ pct_active(s) pct_imagery(s) pct_sensory(s) pct_overlap(s) r(s)])

        if s == 3
            figure(3)
            subplot(2,4,i)
            side = round(sqrt(length(mfr_area_imagery)));
            imagesc(reshape(mfr_area_imagery,side,side));
            colormap(gray);
            title(name{i})

            figure(4)
            subplot(2,4,i)
            side = round(sqrt(length(mfr_area_sensory)));
            imagesc(reshape(mfr_area_sensory,side,side));
            colormap(gray);
            title(name{i})

            figure(5)
            subplot(2,4,i)
            side = round(sqrt(length(mfr_area_sensory)));
            imagesc(reshape(mfr_area_imagery-mfr_area_sensory,side,side));
            %colormap(gray);
            title(name{i})
           
            figure(6)
            subplot(2,4,i)
            side = round(sqrt(length(mfr_area_sensory)));
            imagesc(reshape(active_imagery-active_sensory,side,side));
            colormap(gray);
            title(name{i})
 
        end
        
    end
    disp('---------------------------------------------------------------');
    disp([ mean(pct_active) mean(pct_imagery) mean(pct_sensory) mean(pct_overlap) mean(r)])
    disp(' ' );    
    

end
